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ABSTRACT 


The background response of a gas ionizing proportional X-Ray 
detector is estimated by solving the one dimensional photon 
transport equation for two regions using Monte Carlo techniques . 
The solution was effected using the SSL VAX 780 and the CRAY XMP 
computers at Marshall Space Flight Center. The isotropic photon 
energy spectrum encompassing the range from 1 to 1000 KeV 
incident onto the first region, the shield, is taken so as to 
represent the measured spectrum at an altitude of 3 mb over 
Palastine, Texas. The differential energy spectrum deposited in 
the gas region, xenon, over the range of 0 to 100 KeV is written 
to an output file. In addition, the photon flux emerging from the 
shield region, tin, over the range of 1 to 1000 KeV is also 
tabulated and written to a separate file. Published tabular cross 
sections for photoelectric, elastic and inelastic Compton 
scattering as well as the total absorption coefficient are used. 
Histories of each incident photon as well as secondary photons 
from Compton and photoelectric interactions are followed until 
the photon either is absorbed or exits from the regions under 
consideration. A study is made of the effect of shielding 
thickness upon the energy spectrum deposited in the xenon region 
for this background spectrum incident upon the tin shield. 
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INTRODUCTION 


X-ray sources that are of interest to astrophysical research 
appear very faint when observed from the vicinity of earth. This 
quite low signal intensity puts severe constraints onto the 
detectors, especially imaging types, to be used in planned 
orbiting observatories. In the interest of increasing the signal 
to noise in future detector design by reducing the background, a 
study was undertaken of the background sensitivity of a possible 
type of detector being considered for the AXAF project scheduled 
for the near future. There are advantages to using a position 
sensitive proportional detector and considerable development of 
this type of counter has transpired. As a first step in the study 
of the background to be expected in such a detector from the 
ambient photon flux, a one dimensional geometry calculation of 
tne transport through two regions was made. 

The mathematical description of the transport of so called 
"Hard X-rays" through matter yields integral equations that 
become unwieldy in even the simplest of geometrical arangements . 
Instead of an analytic solution, it is frequently adequate to 
perform a schostacic or Monte Carlo type calculation. In the 
Monte Carlo calculation, individual photon histories are 
catalogued as they transit the various media characterized by its 
cross sections. The result of such a calculation could be either 
the energy spectrum of photons at a surface in the medium or the 
energy deposited due to the transport of photons into a volume of 
the medium. 

Central to a Monte Carlo calculation is a sequence of random 
(or at least pseudo-random) numbers. The sequence used in this 
calculation was generated by the FORTRAN routine RAN(J) on the 
SSL VAX 780 or RANF( ) on the CRAY. In order to test just how 
random this sequence was, 10 million random numbers were 
generated and sorted into 100 intervals evenly distributed from 0 
to 1. The number in each bin, divided by the 10 million, could be 
interpeted as the probability of occurance of a number within 
each interval. Such a distribution is shown in Figure 1 and it 
does indeed appear quite uniform. Figure 2, however, shows the 
deviation from the average of each bin which is a more severe 
test of randomness and, while the deviations are quite small, 
some of the "Pseudo" of the randomness of the sequence may be 
seen . 
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THE INCIDENT X-RAY SPECTRUM 


The incident photon spectrum was taken as: 


dl/dE = 26 (p/3 ) 0 * 3 E -3 ' 97 photons cm”^s“^KeV'^. (1) 


An expression which reasonably describes the energy spectrum 
measured at balloon altitudes over Palastine, Texas (Dean and 
Dipper 1981). The zenith angle of the photons incident upon the 
horizontal plane (the tin sheet) was considered isotropic. At a 
given pressure, the above expression may be considered, when 
properly normalized over the interval [0,1], the probability 
density in the variable E. The prescription for selecting the 
energies at random so that as the number of selections approach 
infinity the energies will be distributed as the above expression 
yields : 


1-r = ( Eq/E ) 0 ’ 97 


( 2 ) 


Where r is a member of a random sequence over the interval [0,1] 
which, in practice, is a pseudo-random sequence generated by a 
computer. Eq is a cutoff energy which is required for 
normalization because the expression does not remain finite as E 
approaches zero. Also since 1-r is itself just a random sequence 
it can be replaced by r. The expression may be solved for E in 
terms of r and then the energies (E) calculated from the 
resulting formula will be distributed according to Equation (1). 
Inverting Equation (2) yields: 


E = (1.0/r) 1 ’ 031 


( 3 ) 


To verify the validity of the Monte Carlo scheme a program was 
written to pick one million energies over the interval 1 to 1000 
KeV using Equation (3) and the RAN(J) random number generator on 
the SSL VAX 780 computer. These energies were sorted into 10 KeV 
intervals to simulate the original spectrum. Equation (1). 
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Since the same distribution can be calculated analytically by 
just integrating Equation (1) over each of the 10 KeV intervals, 
this was done and TABLE 2 shows the comparison of the one million 
history Monte Carlo simulation with the analytic calculation. 


TABLE 1 - COMPARISON OF MONTE CARLO (A MILLION HISTORIES) 
WITH THE ANALYTIC CALCULATION OF 
THE INCIDENT ENERGY SPECTRUM 


Interval 

Monte Carlo 

Analytic 

00.00 -> 

10.99 

23.4613 

24.1833 

10.99 -> 

20.98 

1.8475 

1.2211 

20.98 -> 

30.97 

0.4256 

0.4404 

30.97 -> 

40.96 

0.2210 

0.2278 

40.96 -> 

50.95 

0.1344 

0.1396 

50.95 -> 

60.94 

0.0912 

0.0944 

60.94 -> 

70.93 

0.0668 

0.0681 









930 . 07->940 . 06 

0.0003068 

0.0003648 

940.06->950.05 

0.0003536 

0.0003572 

950 . 05->960 . 04 

0.0002990 

0.0003499 

960.04->970.03 

0.0003510 

0.0003428 

970 . 03->980 . 02 

0.0003198 

0.0003359 

980 . 02->990 . 01 

0.0003276 

0.0003292 

990 . 01->1000 . 0 

0.0002782 

0.0003227 


In the interest of brevity, the entire spectrum from 1 to 
1000 KeV is not presented, only the begining and the end. As may 
be noted by examining TABLE 2, 1 million histories are not 
adequate to represent the desired energy spectrum where the 
derivative of the curve is large (low energies) and where the 
magnitude of the flux is very small (high energies). Figure 3 
shows a plot of the Monte Carlo simulation of the incident 
spectrum at the higher energies which is to be compared with 
Figure 4, the analytic calculation of the incident spectrum, 
integration of Equation ( 1 ) over each interval . Here again the 
inadequacy of a million histories is evident. Consequently, all 
subsequent Monte Carlo calculations made using this 
representation of the incident spectrum were made by compiling 10 
million or more photon histories. 
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THE MEAN FREE PATH 


A photon of given energy while traversing a medium, in 
contrast to a charged particle, does not have a well defined 
range or distance it travels before undergoing any one of a 
number of possible interactions with the medium. The number of 
photons, dN, that have an interaction within an infinitesmal 
distance dx is given by (Jenkins, Gould, and Gedcke, 1981): 


dN = s N dx (4) 


Where s is the total cross section in inverse centimeters, x 
being in centimeters, and is a function of the photon's energy. N 
is the number of photons present at the position x. s is the sum 
of the cross sections for the photoelectric, Compton elastic, and 
Compton inealstic interactions. The energies being considered in 
this work are below the threshold for pair production. 

Equation ( 4 ) may be rearranged as : 


dN/N = s dx 


(5) 


This differential equation has the following solution: 


N/Nq = e -sx 


( 6 ) 


Nq is the flux at x = 0 and the ratio N/Nq is the probability 
that a photon travels the distance x before it undergoes an 
interaction and it is properly normalized for all positive x. 
This probability function may be solved for x using a random 
sequence r as has been done before: 


x = -(1.0/s) logrithm(r) 


(?) 


Where the natural logrithms are to be used. 

A FORTRAN program was written which computed a million 
values of x using Equation (7), with s taken as 1, using the 
function RAN ( J ) to generate a random sequence. These values were 
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sorted into 10 equal intervals over the range of 0 to 1 and 
normalized for comparison with the analytic evaluation of 
Equation (6) over the same intervals. This comparison is shown in 
TABLE 2. 


TABLE 2 - EVALUATION OF THE FREE PATH DISTRIBUTION 
ONE MILLION HISTORIES 


Interval 

Monte Carlo 

Analytic 

0 -> . 1 

.094876 

.0951626 

.1 -> .2 

.086219 

.0861066 

.2 -> .3 

.077935 

.0779126 

.3 -> .4 

.070229 

.0704982 

.4 -> .5 

.063695 

.0637894 

.5 -> .6 

.057681 

.0577191 

.6 -> .7 

.052426 

.0522263 

.7 -> .8 

.047123 

.0472643 

.8 -> .9 

.042736 

.0427513 

.9 -> 1. 

.038735 

.0386902 


It may be noted that the Monte Carlo calculation, compiling 1 
million histories, differs from the analytic calculation by less 
than one percent over the entire range of the table. 
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THE CROSS SECTIONS 


Values of the photoelectric, Compton elastic, Compton 
inelastic, and total cross sections as functions of the photon 
energy, were taken from the literature (Veigele 1973), multiplied 
by the density (Weast 1987), and entered into sequential data 
files for the SSL VAX 780 and the CRAY computers. There is a 
separate file for each element. The file structure, not all 
records , is shown below in TABLE 3 . The fields in the first 
record are number of records in the table, and then the energies 
of the first four absorption edges. Each subsequent record 
contains five fields also: energy, total, photoelectric, Compton 
elastic, and Comption inelastic. The cross section values are in 
inverse centemeters . Logrithmic interpolation in both energy and 
cross section is made for energies between the tabulated values 
and, if the energy at which the cross section is desired falls 
exactly on an absorption edge, the greater value of the cross 
section is used. A separate program was written to check the 
accuracy of the cross section tables' transcription to the files 
from the journals . 


TABLE 3 - STRUCTURE OF XENON CROSS SECTION FILE 


37,4.782,5.104,5.453,34.56 

1.0. 60.593.60.593. .04607. .00002431 
1.145,43.317,43.317, .04531, .00003054 

1.145.47.573.47.573. . 64531,-00003054 

1.5.24.262.24.262. . 04356,-00004682 

2.0. 11.843.11.816. .04056. .00007035 

3.0. 4.332.4.306. .03530. .0001162 

4.0. 2.128.2.098. .03054. .0001560 

4.782.1.372.1.344. . 02754. .0001848 

4.782.40.312.40.06. . 02754. .0001848 

5.0. 3.555.3.530. .02679. .0001923 

5.104.3.380.3.355. . 02654. .0001955 

5.104.4.707.4.682. . 02654. .0001955 

5.453.3.931.3.906. . 02529,-0002065 

5.453.4.732.4.707. . 02529. .0002065 

6.0. 3.656.3.630. .02371. .0002228 

8. 0. 1. 6 80. 1.660.. 01878. .0002754 
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DETECTOR SIMULATION 


The straightforward calculation of the energy deposited in 
the xenon by the incident photon spectrum discussed above after 
it had passed through a layer of tin involves the solution of an 
integral equation wherein the integrand is not a simple function 
(Jenkins, Gould, and Gedcke, 1981). What is frequently done, 
instead of analytically solving the equation representing the 
response of the device, which of course would be generally 
preferred, is to compromise by performing a "Monte Carlo" or 
stochastic simulation of the process (Kalos and Whitlock, 1986), 
(Rubinstein, 1981), and (Hammersley and Handscomb, 1967). The 
first step in the simulation consists of dividing the physical 
process into a number of sequential events whose outcomes are 
governed by known probability distribution functions. The 
outcome: of these events are then determined by properly choosing 

f com these probability distribution functions. The individual 
histories of the objects under study, in this case x-ray photons, 
are then compiled. These historical compilations are the desired 
solution which, in this case is, the differential energy spectrum 
deposited in the xenon layer by a photon spectrum after it first 
passes through a layer of tin. 

For the particular case under study, the sequence of events 
in the history of a photon starts by having its energy selected 
according to Equation (3). Next the direction of the photon, 
polar angle only, is selected so as to reflect the isotropic flux 
onto a horizontal plane. This distribution, after considering the 
solid angle and the projection onto the horizontal plane factors, 
one has: 


P ( x ) = 2x , 


( 8 ) 


where x is the cosine of the polar angle. Equation (8) may be 
integrated to give: 

r = x 2 (9) 

Where r is again a member of a random sequence. Equation (9) may 
be easily solved for the cosine of the polar angle, x, and that 
is what will be required in the program. 

Once an energy and a polar angle for the photon is selected, 
it starts its migration through the tin region by having the tin 
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cross sections written into the active array. The photon history 
in the tin will be discussed first because, although, the 
historical sequence is the same in tin as it is in xenon, 
different historical facts are compiled in the two cases. The 
free path of this photon, the distance to where it interacts with 
the medium, is then computed from Equation (7) using the total 
cross section, s, in this case for tin, and its z coordinate is 
updated. The z coordinate, serves as a tag to determine which 
medium the photon is in at any time, is checked against the 
cumulative thicknesses of the media and, depending upon which 
medium it is in, the type of interaction is selected. The 
selection consists of picking a random number and comparing its 
magnitude with the cumulative fraction of the total cross section 
that the various types of reactions cross sections are at that 
energy. To consider a particular example, let's suppose that the 
total cross section is 100 while, at that same energy, the 
photoelectric, Compton elastic, and Compton inelastic are 80, 15, 
and 5 respectively. Now, if the random number selected is .54321, 
the interaction would be photoelectric; while, if the random 
number selected was .96752, the type of interaction would be 
Compton inelastic. 

If the type of interaction turns out to be Compton elastic, a 
photon with the same energy is started again from this z and the 
polar angle is chosen uniformly over the interval 0 to pi. If, 
however, the interaction was a Compton inelastic interaction, a 
new photon of reduced energy E, is started at that z. The energy 
of the emergent photon from this reaction site is given by 
(Jenkins, Gould, and Gedcke, 1981): 


E-l = E/(1.0 + ( E/511 ) * ( 1 . 0 - cos ( x ) ) ) (10) 


Its polar angle is also chosen uniformly over the interval 0 to 
pi. By choosing the polar angle uniformly here one deviates from 
reality somewhat for the higher energy photons. The Klein-Nashina 
formula (Jenkins, Gould, and Gedcke, 1981) should be used for the 
higher energy photons, which is not uniform in scattering angle 
but peaks in the forward direction. 

Photoelectric interactions give a much wider choice of 
possibilities. If the photon energy is below the L-edge, it is 
considered absorbed in the xenon and the history ends in the tin. 
For energies above the L-edge but below the K-edge, only an 
L-vacancy could be excited in an atom when the photon was 
absorbed photoelectrically . The probability that an L-vacancy was 
formed was calculated from the "jump factors" (Jenkins, Gould, 
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and Gedcke , 1981) (Conde and dos Santos 1985) (McMaster, Kerr, 
Mallet, and Hubbell 1969) at the photon energy by comparing its 
magnitude with a member of the random sequence. This excited ion 
could then either decay radiatively, the emission of a 
fluorescent L photon, or it could decay non-radiatively, i.e. 
Auger processes . Which happened was determined by comparison of a 
random number and the fluorescent yield w K r (Bambynek, Crasemann, 
Fink, Frund, Swift, Price, and Rao 1972). A photon whose energy 
was above the K-edge could create either a K- or L-vacancy in the 
atom that absorbed it . Which it created was determined by 
comparing the magnitude of a member of a random sequence and the 
jump factors for the K and L edges in the same fashion that was 
done in selecting the type of interaction from the cross sections 
as was described above. If the random number's magnitude was 
greater than the combined probability of creating a K- added to 
the probability of creating an L-vacancy, the photon was 
considered absorbed (it really excited a higher level but we are 
not considering them) . If a K-vacancy was formed, the probability 
that it decayed by radiative emission, a fluorescent K photon, or 
non-radiative (Auger process) was taken as the fluorescent yield 
Wg (Bambynek, Crasemann, Fink, Frund, Swift, Price, and Rao 
1972). It was then decided if the K photon emitted was a K-alpha 
or a K-beta using the yield ratios as the relative probability 
(Johnson and White 1985). In case it was a K-alpha, an excited 
ion was left at that position with an L-vacancy which would have 
to be dealt with after the K-alpha 's history terminated. 

In the tin region, a photon's history terminated when it was 
absorbed or escaped from the surface of the tin opposite the 
xenon i.e. was reflected back in the incident direction so as to 
make z less than zero. The energy spectrum of the photons that 
leave the tin into the xenon is compiled and serves as one of the 
output files. 

The history of the photons that enter the xenon are continued 
in the same fashion as they were in the tin except they are 
restarted at the surface of the xenon with the same energy and 
direction they had when they left the tin. Their paths and 
interactions are now governed by the xenon cross sections and the 
absorbed energy spectrum is compiled for them. 
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RESULTS AND CONCLUSION 


The present version of the program provides two output files : 
ENGDEPXE.DAT and ENGEXTIN.DAT. ENGDEPXE.DAT is the energy 
deposited in the xenon region over the range of 0 to 100 KeV in 1 
KeV intervals . ENGEXTIN.DAT is the photon energy spectrum that 
penetrates through the tin shield and enters the xenon region. 
This spectrum covers the energy range of 1 to 1000 KeV in 100 10 
KeV intervals. Figure 5 is a plot of the incident spectrum. 
Equation (3), superimposed over the spectrum exiting a 3/32 inch 
tin shield. It may be seen that, while the tin shield is very 
effective in attenuating the background flux at the lower 
energies, the shield is effectively transparent to the higher 
energy photons . 

Figure 6 is a plot of the spectrum emerging from the tin 
shield and entering the xenon region except the vertical scale is 
expanded to emphasise the lower energy part. It may be noted 
that, except for the fluorescent tin peaks, there are no photons 
below 60 KeV which penetrate the shield. There is only one 
fluorescent K peak shown on the plot because the "resolution" of 
10 KeV per channel does not permit separation of the K-alpha and 
K-beta peaks. There appears, on close inspection, even a very 
small L fluorescent peak in the first channel. 

Figure 7 is a plot of the spectrum of the energy deposited in 
the xenon region (30.0 cm @ 1 atmosphere) for the same shielding 
as was present in Figure 3. Here both the K-alpha and K-beta 
fluorescent peaks of tin are resolved and a slight background is 
present due to higher energy photons that Compton scatter in 
either the shield or the xenon itself. 

Figure 8 is a plot of the background spectrum deposited in 
the xenon with the tin fluorescence "turned off" for several 
thicknesses of the tin shield. It may be noted that for 0.1 cm of 
tin that there is a considerable background. Increasing the 
thickness of tin to 0.238125 (=3/32 inch) reduces this background 
considerably but that additional shielding does not appear to 
reduce the background significantly. 
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APPENDIX I 


FORTRAN PROGRAM TO SIMULATE AN X-RAY DETECTOR 
ON SSL-VAX FILE: [ ELDRIDGE ] PHOTRAN2 . FOR 
ON CRAY FILE: EADS : :HHBE219 . PROGRAMS . FORTRAN( TRAN 2 ) 

PROGRAM PHOTRAN2 

♦COMPUTES THE ENERGY DISTRIBUTION DEPOSITED IN XENON AFTER 

* TRANSPORT THROUGH A TIN PLATE SHIELD. 

* USES SUBROUTINES; INCPHO, XSEC, MFP, XENON 

* INPUT FROM KEYBOARD; TINTH, XENONTH , PRESSURE , N(number of histories) 

REAL *4 EP ( 50 ) , SIGTOT(50) / SIGPHO(50), SIGCCO(50), SIGCIN(50), 

2 EB ( 101 ) , E, SIGT, SIGP, SIGC, SICI, X, INTENSITY ( 101 ) , 

2 XSECT( 2,5,50), EA(101), PHO(lOl) 

INTEGER* 4 J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGCIN, NXC, EB, INTENSITY, 

2 XSECT, TINTH, XENONTH, PRESSURE 

J=9 87654321 ! SEED THE RANDOM NUMBER GENERATOR 

* READ IN FROM THE CONSOLE THE NUMBER OF PHOTONS THAT ARE TO BE DONE 

PRINT * , ' TYPE IN THICKNESS OF THE TIN IN CENTIMETERS ( X . XXX ) ' 

READ ( * , 2 ) TINTH 

PRINT * , ' TYPE IN THICKNESS OF THE XENON IN CENTIMETERS ( X . XXX ) ' 

READ ( * , 2 ) XENONTH 

PRINT * , ' TYPE IN THE XENON PRESSURE IN ATMOSPHERES ( X . XXX ) ' 

READ ( * , 2 ) PRESSURE 

PRINT * , ' TYPE IN THE NUMBER OF PHOTON HISTORIES TO BE DONE (XXX)' ■w' 
READ (*,1) N 

1 FORMAT (18) 

2 FORMAT (F12.6) 

* READ IN THE CROSS SECTION DATA FROM FILES: TINXSECT.DAT, XENXSECT.DAT 

OPEN (UNIT=1 , FILE=' [ELDRIDGEJTINXSECT.DAT' , STATUS='OLD' ) 

OPEN ( UNIT=2 , F ILE= ' [ ELDRIDGE ] XENXSECT . DAT ' , STATUS 3 ' OLD ' ) 

DO 50 L=1 , 2 1 XSECT (L,M,N) L=THE ELEMENT; M=TYPE XSECT; N=ENERGY 

READ (L, 1000 ) XSECT(L, 1,1), XSECT (L, 2 , 1 ) , XSECT(L, 3 , 1 ) , 

XSECT(L, 4,1), 

2 XSECT(L, 5,1) 

NXC = INT ( XSECT ( L , 1 , 1 ) ) 

DO 10 1=2, NXC I READ IN THIS ELEMENT'S CROSS SECTIONS 

READ ( L, 1000 ) (XSECT(L,M, I ) , M=l,5) 

10 CONTINUE 
50 CONTINUE 

1000 FORMAT ( 5F13 . 6 ) 

CLOSE ( UNIT=1 , STATUS 3 ' KEEP ' ) 

CLOSE ( UNIT=2 , STATUS 3 ' KEEP' ) 

* SET UP THE ENERGY BINS FOR THE OUTPUT SPECTRUM 
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EB ( 1 ) = 0.0 
DO 20 I = 2, 101 

EB ( I ) = EB(I-l) + 1.0 
20 CONTINUE 

EA( 1 ) = 0.0 
DO I = 2,101 
EA( I ) = EA(I-l) + 10.0 
ENDDO 

DO 300 I = 1, N 
CALL INCPHO(E) 

IF (E) 100, 100, 400 

* IF ( E .LT. 100.0) GOTO 100 

400 CTHETA = SQRT(RAN( J) ) 

* SET CROSS SECTION TABLES FOR TIN 

NXC = INT ( XSECT (1,1,1)) 

DC IN=2 , NXC 

EP ( IN-1 )=XSECT( 1 , 1 , IN) 

SIGTOT ( IN- 1 ) =XSECT ( 1 , 2 , IN ) 

S IGPHO ( IN- 1 ) =XSECT ( 1 , 3 , IN ) 

SIGCCO( IN-1 )=XSECT( 1,4, IN) 

SIGCIN ( IN-1 ) =XSECT( 1 , 5 , IN ) 

END DO 
NXC=NXC-1 

500 CALL XSEC ( E , SIGT, SIGP, SIGC, SIGI) 

CALL MFP ( SIGT , PATHLEN) 

Z = Z + PATHLEN*CTHETA 

IF ( Z .GT. TINTH .AND. Z .LE. (TINTH+XENONTH) ) THEN 
DO III = 2,101 

I F ( E . GE . EA ( I I I - 1 ) . AND . E . LT . EA (III) ) PHO ( III-l ) =PHO ( I I 1-1 ) + 1 . 0 
ENDDO 

CALL XENON ( E , Z, CTHETA} 

GOTO 1100 

ELSE IF (Z .LE. 0.0 .OR. Z .GT. (TINTH+XENONTH)) THEN 
GOTO 1100 
ELSE 

CONTINUE 
END IF 

* PHOTON IS STILL INSIDE THE PLATE; SELECT THE TYPE OF INTERACTION 

* SET CROSS SECTION TABLES FOR TIN 

NXC = INT ( XSECT (1,1,1)) 

DO IN=2 , NXC 

EP ( IN- 1 ) =XSECT ( 1 , 1 , IN ) 

SIGTOT ( IN-1 ) =XSECT (1,2, IN ) 

S IGPHO ( IN- 1 ) =XSECT ( 1 , 3 , IN ) 

S IGCCO ( IN- 1 ) =XSECT ( 1 , 4 , IN ) 

SIGCIN ( IN- 1)=XSECT( 1,5, IN) 
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200 


END DO 
NXC=NXC-1 

X = RAN ( J ) 

Y = SIGP/SIGT 
W = SIGC/SIGT 

V = SIGI/SIGT 

IF (0.0 .LE. X .AND. X .LT. Y) GOTO 600 
IF (Y .LE. X .AND. X .LT. (Y+W) ) GOTO 700 

* PICK ANGLE (NOT — KLEIN-NASHINA) ! COMPTON INELASTIC 

800 THETA = 3 . 14159 27*RAN(J) 

CTHETA= COS (THETA) 

* COMPUTE NEW ENERGY DEPENDS ON ANGLE 

E « E/( 1.0+ (E/511 )*(1. 0-CTHETA) ) 

GOTO 500 

* PHOTOELECTRIC INTERACTION 

600 IF ( E .LT. 3.929) GOTO 1100 

IF ( E .LT. 29.2) GOTO 704 

* BOTH K AND L VACANCIES ARE ENERGETICALLY POSSIBLE 

X = RAN ( J ) 

IF ( X .LT. 0.84556) GOTO 701 

703 IF ( X .LT. 0.949529) GOTO 702 

* GOT NEITHER K NOR L VACANCY THIS TIME; START NEW PHOTON 

GOTO 1100 

* GOT AN L VACANCY; DID IT DECAY BY RADIATIVE TRANSITION? 
702 X = RAN(J) 

IF ( X .GT. 0.064) GOTO 1100 

* YES IT WAS A RADIATIVE TRANSITION; ENERGY OF 3.44 KEV 

E = 3.44 
GOTO 700 

* GOT A K VACANCY; DID IT DECAY BY RADIATIVE TRANSITION? 
701 X = RAN ( J ) 

IF ( X .GT. 0.859) GOTO 1100 

* YES , IT WAS A RADIATIVE TRANSITION; K-ALPHA OR K-BETA?? 

X=RAN( J ) 

IF ( X .LE. 0.78) THEN 
E = 25.2 ! IT WAS A K-ALPHA 

KL “ KL + 1 ! CREATES AN L VACANCY 

ELSE 

E = 28.5 
END IF 
GOTO 700 

* ONLY AN L VACANCY WAS ENERGETICALLY POSSIBLE 

704 X = RAN ( J ) 

IF ( X .GT. 0.6732) GOTO 1100 

* GOT AN L VACANCY; DID IT DECAY BY RADIATIVE TRANSITION? 

GOTO 702 
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* ALSO COME TO HERE FOR ELASTIC COMPTON EVENTS 
700 THETA = 3 . 1415927*RAN( J) 

CTHETA = COS (THETA) 

GOTO 500 

* PHOTON TRAVERSED THE PLATE; START NEXT PHOTON AT THE SURFACE 

* FIRST, HOWEVER; ARE THERE ANY L VACANCIES LEFT?? 

1100 IF(KL .NE. 0) THEN 

KL = KL - 1 I FINISH OFF THE L-VACANCIES CREATED EARLIER 
GOTO 702 
END IF 

* GO BACK AND START ANOTHER PHOTON 
100 Z = 0.0 

300 CONTINUE 

* UNLESS; OF COURSE, FNINISHED; WRITE THE SPECTRUM TO A FILE 

OPEN (UNIT=1, FILE=' [ELDRIDGE] ENGDEPXE.DAT ' , STATUS= 'OLD') 
OPEN (UNIT=2, FILE=' [ELDRIDGEJENGEXTIN.DAT' , STATUS= 'OLD') 

DO 40 L = 1, 100 

* NORMALIZE INTENSITY; FACTOR OF 26.0 FROM ENERGY DISTRIBUTION 

INTENSITY (L) = 26 . 0*INTENSITY(L) /FLOAT(N) 

PHO(L) = 26.0* PHO ( L ) / FLOAT ( N ) 

WRITE (1,900) EB(L), INTENSITY(L) 

WRITE (2,900) EA(L), PHO(L) 

900 FORMAT (F8.2, F20.10) 

40 CONTINUE 

CLOSE ( UNIT=1 , STATUS='KEEP' ) 

CLOSE ( UNIT=2 , STATUS =' KEEP ' ) 

* THEN, AT LAST, STOP 1 ! ! ! 1 

STOP 

END 


SUBROUTINE INCPHO(E) 

* SELECTS ENERGY OF PHOTON FROM BACKGROUND DIST. 

* RETURNS A VALID ENERGY FOR A PHOTON IF SUCCESSFUL; HISTORY BEGINS 

* RETURNS 0 FOR PHOTON ENERGY IF NOT SUCCESSFUL; END OF HISTORY 

REAL *4 EP ( 50 ) , SIGTOT(50), SIGPHO(50), SIGCCO(50), SIGCIN(50), 
2 E, SIGT, SIGP, SIGC, SICI, X 

INTEGER*4 J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGCIN, NXC 

* RANDOM ENERGY UNIFORM 0 TO 1000 KEV 

X=RAN ( J ) 

IF ( X .EQ. 0.0) X=RAN ( J ) 

E=1 . 0/X**l . 031 

* DISREGARD PHOTONS WITH: 1 KEV > ENERGY > 1000 KEV 

IF ( E .LT. 1.0 .OR. E .GT. 1000.0) THEN 
E = 0.0 
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RETURN 

* SUCCESS; SO ACCEPT THIS ENERGY; DO PHOTON HISTORY 
ELSE 

RETURN 
END IF 
END 


SUBROUTINE MFP(SIGT, PATHLEN) 

* PICKS THE PHOTON PATH LENGTH FROM THE PROPER DISTRIBUTION 

* WHEN FURNISHED THE TOTAL CROSS SECTION 

REALM EP ( 50 ) , SIGTOT(50), SIGPHO(50), SIGCCO(50), SIGCIN(50), 
2 E, SIGT, SIGP, SIGC, SICI, X 

INTEGER* 4 J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGCIN, NXC 
X = RAN(J) 

PATHLEN = - ( 1 . 0 /S IGT ) *ALOG ( X ) 

RETURN 

END 

* 


* 


SUBROUTINE XSEC(E, SIGT, SIGP, SIGC, SIGI) 

* INTERPLOATES TABLES OF CROSS SECTIONS VERSUS ENERGY 

* WHEN FURNISHED WITH AN ENERGY FROM 1 TO 1000 KeV 

* TOTAL, PHOTOELECTRIC, COMPTON ELASTIC, COMPTON INELASTIC ALL RETURNE 

REALM EP ( 50 ) , SIGTOT(50), SIGPHO(50), SIGCCO(50), SIGCIN(50), 

2 E, SIGT, SIGP, SIGC, SICI, X 

INTEGERM J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGCIN, NXC 
IF (E .LT. EP ( 1 ) .OR. E .GT. EP(NXC)) THEN 
*IS ENERGY OUT OF RANGE ? 

GOTO 500 

* RETURN TO CALLING ROUTINE 

ELSE 

* INTERPOLATE TABLES TO GET X-SECTIONS FOR THAT ENERGY 

* CASE THE ENERGY IS EXACTLY THE FIRST POINT IN THE TABLE 

IF ( E .EQ. EP ( 1 ) ) THEN 
SIGT=SIGTOT( 1 ) 

SIGP=SIGPHO( 1 ) 

SIGC=SIGCCO( 1 ) 

SIGI=SIGCIN ( 1 ) 

GOTO 500 
END IF 

* CASE IF THE ENERGY FALLS EXACTLY ON ANY POINT IN THE TABLE 
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DO 10 K=2 ,NXC 

IF (E .EQ. EP(K) .AND. EP(K) .EQ. EP(K+1)) THEN 
SIGT=SIGTOT(K+l ) l AT A POINT; AN EDGE 
S IGP=S IGPHO ( K+l ) 

SIGC=SIGCCO(K+l ) 

SIGI=SIGCIN(K+1) 

GOTO 500 

ELSE IF (E .EQ. EP(K)) THEN 1 AT A POINT IN TABLE 
SIGT=SIGTOT(K) 

S IGP=S IGPHO (K) 

SIGC=SIGCCO(K) 

SIGI=SIGCIN(K) 

GOTO 500 
ELSE 

OTHERWISE INTERPOLATE TO GET CROSS-SECTIONS AT THAT ENERGY 
IF (E .GT. EP(K-l) .AND. E .LT. EP(K)) THEN 
X=LOG( E/EP (K-l ) )/LOG(EP(K)/EP(K-l) ) 

SIGT=EXP (LOG( SIGTOT(K) / SIGTOT (K-l ) ) *X+LOG( SIGTOT( K-l ) ) ) 
SIGP=EXP (LOG( SIGPHO(K) / SIGPHO(K-l ) ) *X+LOG( S IGPHO (K-l ) ) ) 
S IGC=EXP ( LOG ( S IGCCO ( K ) /SIGCCO ( K- 1 ) ) *X+LOG ( SIGCCO (K-l ) ) ) 
SIG1'=EXP (LOG( SIGCIN(K) /SIGCIN(K-1 ) ) *X+LOG( SIGCIN(K-1 ) ) ) 
END IF 
END IF 

10 CONTINUE 

END IF 
500 RETURN 
END 


SUBROUTINE XENON ( E, Z, CTHETA) 

REAL*4 EP ( 50 ) , SIGTOT(50), SIGPHO(50) / 810000(50), SIGCIN(50), 

2 EB(101), E, SIGT, SIGP, SIGC, SICI, X, INTENSITY (101), 

2 XSECT( 2,5,50), ETOT 

INTEGER*4 J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGCIN, NXC, EB, INTENSITY, 
2 XSECT, TINTH, XENONTH, PRESSURE 

SET CROSS SECTION TABLES FOR XENON 
NXC = INT (XSECT (2,1,1 ) ) 

DO IN=2 , NXC 

EP ( IN-1 ) =XSECT ( 2 , 1 , IN) 

SIGTOT ( IN-1 ) =XSECT( 2,2, IN ) *PRESSURE 
S IGPHO ( IN- 1 ) =XSECT ( 2 , 3 , IN ) *PRESSURE 
SIGCCO ( IN- 1 ) =XSECT ( 2 , 4 , IN ) *PRESSURE 
S IGCIN ( IN- 1 ) =XSECT ( 2 , 5 , IN ) *PRESSURE 
END DO 



NXC=NXC-1 
Z = TINTH 
ETOT = 0.0 

500 CALL XSEC ( E, SIGT, SIGP, SIGC, SIGI) 

CALL MFP(SIGT, PATHLEN) 

Z = Z + PATHLEN*CTHETA 

IF (Z .LE. TINTH .OR. Z .GT. (TINTH+XENONTH) ) GOTO 1100 
200 X = RAN(J) ! IN XENON; PICK TYPE INTERACTION 

Y = SIGP/SIGT 
W = SIGC/SIGT 

V = SIGI/SIGT 

IF (0.0 .LE. X .AND. X .LT. Y) GOTO 600 l PHOTOELECTRIC 
IF (Y .LE. X .AND. X .LT. (Y+W)) GOTO 700 ICOMPTON ELASTIC 
800 THETA = 3 . 1415927*RAN( J) ICOMPTON INELASTIC 

CTHETA= COS (THETA) l ISOTROPIC (NOT VERY GOOD) 

EG - E/ (1.0+( E/511 )*( 1. 0-CTHETA) ) I NEW ENERGY DEPENDS ON ANGLE 
EE = E-EG I THE ELECTRON ENERGY (ABSORBED) 

ETOT = ETOT + EE I TALLY THE DEPOSITED ENERGY 

E = EG I THE NEW ENERGY OF THE PHOTON 

GOTO 500 l START A PHOTON FROM THIS POSITION AT NEW ENERGY 

* PHOTOELECTRIC INTERACTION 

600 IF ( E .LT. 4.781) THEN l ENERGY TOO LOW TO EXCITE L-VACANCY 

ETOT = ETOT + E I THIS PHOTON ABSORBED; HIGHER LEVEL EXCITATION 
GOTO 1100 l ANY MORE EXCITED ATOMS ??? 

END IF 

IF(E .LT. 34.582) GOTO 704 IL-VACANCY ENERGETICALLY POSSIBLE 

* BOTH K AND L VACANCIES ARE ENERGETICALLY POSSIBLE 

X = RAN(J) 

IF(X .LT. 0.8355) GOTO 701 1 GOT A K-VACANCY 

703 IF(X .LT. 0.949529) GOTO 702 1 GOT AN L-VACANCY 

* GOT NEITHER K NOR L VACANCY THIS TIME; START NEW PHOTON 

* THE PHOTON WAS COMPLETELY ABSORBED 

ETOT = ETOT + E 
GOTO 1100 

702 X = RAN(J) 1 GOT AN L-VACANCY; DID IT DECAY BY RADIATIVE TRANSITIO 
IF ( X .GT. 0.1) THEN 

ETOT = ETOT + E 1 NO; COMPLETELY ABSORBED 
GOTO 1100 
END IF 

ETOT = ETOT + E - 4.781 I YES; RADIATIVE 

E = 4.781 I ENERGY EMITTED PHOTON 

GOTO 700 1 START 4.781 KEV PHOTON OUT FROM HERE 

* GOT A K VACANCY; DID IT DECAY BY RADIATIVE TRANSITION? 

701 X = RAN(J) 

IF ( X .GT. 0.88) THEN 
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ETOT = ETOT + E 1 NOT RADIATIVE; COMPLETELY ABSORBED HERE 
GOTO 1100 I ANY MORE EXCITED IONS??? 

ELSE 

* YES, IT WAS A RADIATIVE TRANSITION; WAS IT A K-ALPHA OR K-BETA?? 

X = RAN(J) 

IF ( X .LE. 0.7667) THEN 

ETOT = ETOT + E - 34.582 ! THIS MUCH ENERGY WAS DEPOSITED 

E = 29.801 1 WHEN A K-ALPHA PHOTON WAS EMITTED 

KL = KL + 1 I ALSO AN L-VACANCY PRODUCED 

ELSE 

ETOT = ETOT +E - 33.644 ! THIS MUCH ENERGY WAS DEPOSITED 

E = 33.644 l WHEN A K-BETA PHOTON WAS EMITTED 

END IF 
END IF 
GOTO 700 

704 X - RAN ( J ) ! ONLY AN L-VACANCY WAS ENERGETICALLY POSSIBLE 

IF ( X .GT. 0.6527) THEN 

ETOT -•= ETOT + E ! NO L-VACANCY, HIGHER LEVELS?? 

GOTO 1100 ! ANY MORE EXCITED IONS?? 

END IF 

GOTO 702 ! L-VACANCY WAS PRODUCED; BUT WAS IT RADIATIVE?? 

* ALSO COME TO HERE FOR ELASTIC COMPTON EVENTS 
700 THETA = 3 . 1415927*RAN( J) 

CTHETA = COS (THETA) 

GOTO 500 

* PHOTON TRAVERSED THE XENON REGION; START NEXT PHOTON 

* BUT FIRST; ARE THERE ANY L VACANCIES LEFT IN THE XENON?? 

1100 IF(KL .NE. 0) THEN 

E = 4.781 
KL = KL - 1 
GOTO 702 
END IF 

CALL BINIT (ETOT) 

RETURN 

END 


SUBROUTINE BINIT(E) 

REALM EP ( 50 ) , SIGTOT(50), SIGPHO(50), SIGCCO(50), SIGCIN(50), 

2 EB(101), E, SIGT, SIGP, SIGC, SICI, X, INTENSITY ( 101 ) , 

2 XSECT( 2,5,50), ZL(500) 

INTEGERM J 

COMMON J, EP, SIGTOT, SIGPHO, SIGCCO, SIGC IN, NXC, EB, INTENSITY 
* SELECT THE PROPER BIN FOR THIS PHOTON 
DO 30 K = 2, 101 
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IF (E .GT. EB(K-l) .AND. E .LE. EB(K) ) 
INTENSITY ( K-l ) =INTENSITY ( K-l ) +1 . 0 
ELSE 

CONTINUE 
END IF 
30 CONTINUE 
RETURN 
END 


THEN 
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